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Abstract. We consider collisional models for granular particles and analyze the 
conditions under which the restitution coefficient might be a constant. We show 
that these conditions are not consistent with known coilision laws. From the gener- 
alization of the Hertz contact law for viscoelastic particles we obtain the coefficient 
of normal restitution e as a function of the normal component of the impact veloc- 
ity Vimp. Using e(tiiDip) we describe the time evolution of temperature and of the 
velocity distribution function of a granular gas in the homogeneous cooling regime, 
where the particles collide according to the viscoelastic law. We show that for the 
studied systems the simple scaling hypothesis for the velocity distribution function 
is violated, i.e. that its evolution is not determined only by the time dependence 
of the thermal velocity. We observe, that the deviation from the Maxwellian dis- 
tribution, which we quantify by the value of the second coefficient of the Sonine 
polynomial expansion of the velocity distribution function, does not depend on time 
monotonously. At first stage of the evolution it increases on the mean-collision time- 
scale up to a maximum value and then decays to zero at the second stage, on the 
time scale corresponding to the evolution of the granular gas temperature. For gran- 
ular gas in the homogeneous cooling regime we also evaluate the time-dependent 
self-diffusion coefficient of granular particles. We analyze the time dependence of the 
mean-square displacement and discuss its impact on clustering. Finally, we discuss 
the problem of the relevant internal time for the systems of interest. 



1 Introduction 

Granular gases, i.e. systems of inelastically colliding particles, are widely 
spread in nature. They may be exemplified by industrial dust or interter- 
restrial dust; the behavior of matter in planetary rings is also described in 
terms of the granular gas dynamics. As compared with common molecular 
gases, the steady removal of kinetic energy in these systems due to dissi- 
pative collisions causes a variety of nonequilibrium phenomena, which have 
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been very intensively studied (e.g. p]— |Tl[]). In most of these studies the coeffi- 
cient of restitution, which characterizes the energy lost in the collisions, was 
assumed to be constant. This approximation, although providing a consider- 
able simplification, and allowing to understand the main effects in granular 
gas dynamics, is not always justified (see also the paper by Thornton in 
this book p^|). Moreover, sometimes it occurs to be too crude to describe 
even qualitatively the features of granular gases. Here we discuss the prop- 
erties of granular gases consisting of viscoelastically colliding particles which 
implies an impact- velocity dependent restitution coefficient. The results are 
compared with results for gases consisting of particles which interact via a 
constant restitution coefficient and we see that the natural assumption of 
viscoelasticity leads to qualitative modifications of the gas properties. 
The following problems will be addressed: 

• Why does the restitution coefficient e depend on the impact velocity Wi mp ? 

• How does it depend on the impact velocity? 

• What are the consequences of the dependence of e on Vi mp on the collective 
behavior of particles in granular gases? In particular how does e = e (^i mp ) 
influence: 

• the evolution of temperature with time? 

• the evolution of the velocity distribution function with time? 

• the self-diffusion in granular gases? 

In what follows we will show that the dependence of the restitution coef- 
ficient on the impact velocity is a very basic property of dissipative particle 
collisions, whereas the assumption of a constant restitution coefficient for 
the collision of three-dimensional spheres may lead to a physically incorrect 
dependence of the dissipative force on the compression rate of the colliding 
particles. From the Hertz collision law and the general relation between the 
elastic and dissipative forces we deduce the dependence of the restitution 
coefficient on the impact velocity, which follows purely from scaling consider- 
ations. We also give the corresponding relation obtained from rigorous theory. 
Using the dependence e(i>i mp ) we derive the time dependence of the temper- 
ature, the time-evolution of the velocity distribution function and describe 
self-diffusion in granular gases in the homogeneous cooling regime. 

2 Dependence of the restitution coefficient on the 
impact velocity 

The collision of two particles may be characterized by the compression £ and 
by the compression rate £, as shown on Fig. [l]. The compression gives rise to 
the elastic force F c i(£), while the dissipative force -Fdiss(£,£) appears due to 
the compression rate. 
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Fig. 1. Sketches of two colliding spheres. The compression £ is equal to 2R— \r i— r-2 \, 
with ri/2 being the particles positions. The compression rate is £ = vi — «2- For 
simplicity the head on collision of identical spheres is shown. 



If the compression and the compression rate are not very large, one can 
assume the dependence of the elastic and dissipative force on £ and £ 

Fel(0 ~ £ Q (1) 

l^-O- ^C . (2) 

The dimension analysis yields the following functional form for the depen- 
dence of the restitution coefficient on the impact velocity |l3| : 



<v imp ) = e (v^ +P j (3) 

Therefore, the condition for a constant restitution coefficient imposes the 
relation between the exponents a, (3 and 7 



2( 7 -a)+/3 (1 + a) = 0. (4) 

For compressions which do not exceed the plasticity threshold, the parti- 
cle's material behaves as a viscoelastic medium. Then it may be generally 
shown [|]~6| — that the relation 



F diss =AZ—F el (Q (5) 

between the elastic and dissipative force holds, independently on the shape 
of the bodies in contact, provided three conditions are met JlSj ] : 

(i) The elastic components of the stress tensor cr^ depend linearly on the 
components of the deformation tensor itjfe ]20| ] . 

(ii) The dissipative components of the stress tensor qffi depend linearly on 
the components of the deformation rate tensor Ujfc p0| . 
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(iii) The conditions of quasistatic motion are provided, i.e. £ <C c, t v ; s <C t c 
[16, ^| (here c is the speed of sound in the material of particles and t v - 1s 



is the relaxation time of viscous processes in its bulk) . 
The constant A in Eq. (g) reads 



A= - 



1 (3?7 2 - 



3 (3^2 + 27/!) 



y^ 2 



(6) 



where y and i/ are respectively the Young modulus and the Poisson ratio 
of the particle material and the viscous constants rji, r\ 2 relate (linearly) the 
dissipative stress tensor a^ iss to the deformation rate tensor un- |l6|, fi7| , p0| ]. 
From Eq. (|) follows that 

(3 = 1 7 = a - 1 . (7) 

Consider now a collision of three-dimensional spherical particles of radii 
i?i and i?2- The Hertz contact contact law gives for the elastic force |^l[ 

p cl = P e /2 , P = 2Y ^ vw*, (8) 

3(1 — v z ) 

where R cS = R 1 R 2 /(Ri + R 2 )- With the set of exponents, a = 3/2, (3 = 1 
and 7 = 1/2, which generally follows from the basic laws of the viscoelastic 
collision, the condition for the constant restitution coefficient, Eq. (Q), is 
obviously not satisfied. For spherical particles the restitution coefficient could 
be constant only for 7 = 1/4; this, however, is not consistent with the collision 
laws. Instead one obtains the functional dependence 

^ 5 P ) • 0) 

Note that this conclusion comes from the general analysis of viscoelastic 
collisions with no other assumptions needed. Therefore, the dependence of the 
restitution coefficient on the impact velocity, Eq. (|9|), is a natural property, 
provided the assumption on viscoelasticity holds true which is the case in 



a wide range of impact velocities (see discussion on page pL02| ) . We want to 
mention that the functional dependence Eq. @ was already given in [p2| 
using heuristic arguments. 

Rigorous calculations |^3| yield for the dependence of the restitution co- 
efficient on the impact velocity: 

e = 1 - C lA ^v^ + C 2 A^^vfH T • • • (10) 

with 
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where m = m^/fmi + 1^2) ( m i/2 are the masses of the colliding parti- 
cles) . Numerical values for the constants C\ and C 2 obtained in [^3| may be 
also written in a more convenient form fllal: 



ft - - O.-i0l. (12) 

Although the next-order coefficients of the above expansion C3 = —0.483582, 
C4 = 0.285279, are now available ]13|| , we assume that the dissipative constant 
A is small enough to ignore these high-order terms. (For large A a very 
accurate Pade approximation for e(i>i mp ) has been proposed recently p3| ). 



3 Time-evolution of temperature and of the velocity 
distribution function 



We consider a granular gas composed of N identical particles confined in 
a volume fl. The particles are assumed to be smooth spheres, so that the 
collision properties are determined by the normal component of the relative 
motion only. The gas is supposed to be dilute enough so that one can assume 
binary collisions (i.e. neglect multiple collisions) and ignore the collision du- 
ration as compared with the mean free time in between successive collisions. 
We assume that the initial velocities of the particles (more precisely the tem- 
perature, which we define below) are not very large to assure viscoelastic 
properties of the collisions, i.e. to avoid plastic deformations and fragmen- 
tation. The final velocities are assumed not to be very small which allows 
to neglect surface forces as adhesion and others. Under these restrictions 
one can apply the viscoelastic collision model. Furthermore, we assume that 
dissipation is not large, so that the second-order expansion ( |To| ) for e(t>i mp ) 
describes the collisions accurately. We analyze the granular gas in the regime 
of homogeneous cooling, i.e. in the pre-clustering regime, when the gas is 
homogeneously distributed in space. 

The impact velocity t>i mp of colliding smooth spheres, which determines 
the value of the restitution coefficient according to Eq. (|l(J), is given by the 
normal component of the relative velocity 

«im P = ■ e\ with v 12 = vi-v 2 . (13) 

The unit vector e = t"i2/|t"i2| gives the direction of the intercenter vector 
Ti2 = ri — v<x at the instant of the collision. 

The evolution of the granular gas proceeds by elementary collision events 
in which the pre-collisional velocities of colliding particles vi, v 2 , are con- 
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verted into after-collisional ones, v\, v 2 , according to the rules 

v*=vi--[l + e(|«i2 • e|)] («i2 • e)e 

I (14) 
w 2 = ^2 + ^ I 1 + e (l v i2 ' e|)] (»i2 • e)e , 

where e depends on the impact velocity v- lmp = \v\ 2 ■ e|. Due to the direct 
collision ( |l4| ) the population in the velocity phase-space near the points v\, 
v 2 decreases, while near the points v\, v 2 it increases. The decrease of the 
population near V\, v 2 caused by the direct collision, is (partly) counter- 
balanced by its increase in the inverse collision, where the after-collisional 
velocities are V\, v 2 with the pre-collisional ones v\* , v 2 *. The rules for the 
inverse collision read 

vi = vr-hl + e(\v%-e\)} Kl-e)e 

? (15) 
V2=v* 2 * + -[\ + e{\vH-e\)] K|-e)e. 

Note that in contrast to the case of e = const, the restitution coefficients in 
the inverse and in the direct collisions are different. 

The Enskog-Boltzmann equation describes the evolution of the popula- 
tion of particles in the phase space on the mean-field level. The evolution is 
characterized by the distribution function f(r,v,t), which for the force-free 
case does not depend on r and obeys the equation 0,124) 



= ff2(e0<7 2 J dv 2 J de0(— v 12 ■ e)\v 12 



d_ 

dt" 

x { X f(vT ,t)f(v? ,t) - /(»i,i)/(f7 2 , t)} = g 2 (a)I(f, f) (16) 



where a is the diameter of the particles. The contact value of the pair distri- 
bution function ]2j| 

g 2 (a) = (2- V )/2(l- V ) 3 , (17) 



accounts on the mean-field level for the increasing frequency of collisions due 
to excluded volume effects with r\ = g nna 3 being the volume fraction. 

The first term in the curled brackets in the right-hand side of Eq. ( [l6| ) 
refers to the "gain" term for the population in the phase-space near the 
point Vi, while the second one is the "loss" term. The Heaviside function 
&(—Vi2-e) discriminates approaching particles (which do collide) from sep- 
arating particles (which do not collide), and \v%2 ■ e| gives the length of the 
collision cylinder. Integration in Eq. ( |l6| ) is performed over all velocities v 2 
and interparticle vectors e in the direct collision. Equation ( |l6| ) accounts also 
for the inverse collisions via the factor x, which appears due to the Jacobian 
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of the transformation ■?;**, — > 1)1,1)2, and due to the difference between 
the lengths of the collision cylinders of the direct and the inverse collision: 

V(vi,V2) \vi2-e\ 

For constant restitution coefficient the factor x is a constant 

1 , . 

X = ~2 = const , (19) 

while for e = e (vi mp ), as given in Eq. jTc|), it reads Q 

X = 1 + yC X A« 2 / 5 |t, la • e\^o + |^ s 4/5| Wu . e |2/5 + . . . (20) 



From Eq. ( |20| ) it follows that x = x(l v i2 • Since the average velocity in 
granular gases changes with time, such a dependence of x means, as we will 
show below, that x and, therefore, the velocity distribution function itself 
depend explicitly on time. The time dependence of x changes drastically the 
properties of the collision integral and destroys the simple scaling form of 
the velocity distribution function, which holds for the case of the constant 
restitution coefficient (e.g. ^,||])- 

Nevertheless, some important properties of the collision integral are pre- 
served. Namely, it may be shown that the relation 

j t m)) = J dvuKv^ffat) = 9 2{°) I dvi^Itf, /) = (21) 

g2 ^ )fj2 J dv x dv 2 J deO(-v 1 2-e)\v 1 2-e\f(v 1 ,t)f(v2,t)A[4 ! (v 1 )+yj(v 2 )} 
holds true, where 

(^(t)) = J dwl>(v)f(v,t) (22) 

is the average of some function ip(v), and 

A^{ Vi ) = [^(<) - ^{Vi)\ (23) 

denotes the change of ip( v i) m a direct collision. 

Now we introduce the temperature of the three-dimensional granular gas, 

\nT{t) = J dv^f(v,t), (24) 

where n is the number density of granular particles (n — N/fi), and the 
characteristic velocity Vq (t) is related to temperature via 

T{t) = ±mvl{t) . (25) 
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First we try the scaling ansatz 

^ = W) Rc) (26) 

where c = v/vo(t) and following (J,0| assume that deviations from the 
Maxwellian distribution are not large, so that /(c) may be expanded into 
a convergent series with the leading term being the Maxwellian distribution 
cj){c) = 7T -3 / 2 exp(— c 2 ). It is convenient to use the Sonine polynomials expan- 
sion |,0 

/>) = <Mc) jl + f> P Mc 2 )j • (27) 
These polynomials are orthogonal, i.e. 

dc^{c)S p (c 2 )S pl (c 2 ) = 8 m ,N p , (28) 



where S pp i is the Kronecker delta and M p is the normalization constant. The 
first few polynomials read 



S (Z) = 1 
S 2 (x) 



2 (29) 
x 2 5x 15 



2 2 8 

Writing the Enskog-Boltzmann equation in terms of the scaling variable C\ , 
one observes that the factor \ may not be expressed only in terms of the 
scaling variable, but it depends also on the characteristic velocity vo(t), and 
thus depends on time. Therefore, the collision integral also occurs to be time- 
dependent. As a result, it is not possible to reduce the Enskog-Boltzmann 
equation to a pair of equations, one for the time evolution of the temperature 
and another for the time-independent scaling function, whereas for e =const. 
the Boltzmann-Enskog equation is separable, e.g. @,|5[0|. Formally adopting 
the approach of Refs. j|,|7j for e =const., one would obtain time-dependent co- 
efficients a p of the Sonine polynomials expansion. This means that the simple 
scaling ansatz (^) is violated for the case of the impact-velocity dependent 
restitution coefficient. 

Thus, it seems natural to write the distribution function in the following 
general form 



/(„,*) = _/(c,i) (30) 
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with 

/(c) = 0( c )|l + fj ap (t)5 ?) ( C 2 )| (31) 

and find then equations for the time- dependent coefficients a p (t). Substituting 
( |30|) into the Boltzmann equation (16) we obtain 



M2 
3 



(32) 



with 

B = B(t) = v (t)g 2 (<j)<T 2 n. (33) 
We define the dimensionless collision integral: 

/(/,/> 

J dc 2 J de0(-c 12 ■ e )|ci2 • e| {x/^*,*)/^*,*) - /(<*,*)/(<*,*)} (34) 

with the reduced factor \ 

X = 1 + yCi«5'|ci2 • e| 1/5 + ||c?tf' a | Cia • e| 2 / 5 + • • • (35) 

which depends now on time via a quantity 

<5'(t) = Ak 2/5 [2T(i)] 1/10 = 6 [2T{t)/T } 1/w . (36) 

Here S = Aft 2 / 5 ^] 1 / 10 , To is the initial temperature, and for simplicity we as- 
sume the unit mass, m = 1. We also define the moments of the dimensionless 
collision integral 

lh = ~j dcK?j(fj) , (37) 

so that the second moment describes the rate of the temperature change: 

§--!«■«. (38) 

Equation (Q) follows from the definitions of the temperature and of the 
moment \x 2 . Note that these moments depend on time, in contrast to the 
case of the constant restitution coefficient, where these moments are time- 
independent H. 
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Multiplying both sides of Eq. ( |3^ ) with and integrating over dci, we 
obtain 

oo 

y-p(c p ) - B~ 1 '^a k v kp = n p (39) 
fe=i 

where integration by parts has been performed and we define 

v kv = I ' cp{c)J ] S k (c 2 )dc (40) 

C Pf(c,t)dc. (41) 



The calculation of v kp is straightforward; the first few of these read: ^ 22 = 0, 
1/24 = X' The odd moments <^c 2n+1 ) vanish, while the even ones (c 2n ) may 
be expressed in terms of a k with < k < n, namely, (c 2 ^ = | — |ai. 
On the other hand, from the definition of temperature and of the thermal 
velocity in Eqs. ( p5| ) and (24) follows that (c 2 ) — | and thus, a\ = 0. Similar 



considerations yield (c 4 ) — (1 + 02). The moments /i p may be expressed 
in terms of coefficients 0,2, 0,3, ■■ ■ too; therefore, the system Eq. ( |39| ) is an 
infinite (but closed) set of equations for these coefficients. 

It is not possible to get a general solution of the problem. However, since 
the dissipative parameter S is supposed to be small, the deviations from the 
Maxwellian distribution are presumably small too. Thus, we assume that 
one can neglect all high-order terms with p > 2 in the expansion (plT) . Then 



Eq. (^9|) is an equation for the coefficient a 2 . For p = 2 Eq. (|39|) converts 
into an identity, since (c 2 ) = |, a\ = 0, v 2 2 = and v 2 a = x- F° r P = 4 we 
obtain 

a 2 - - Bfi 2 (1 + a 2 ) + —Bin = . (42) 
6 15 

In Eq. (^) B depends on time as 

B(t) = (STr)- 1 / 2 ^)- 1 ^)/^] 172 , (43) 

where r c (0) is related to the initial mean-collision time, 

r^O)- 1 = An 1 l 2 9 2{a)a 2 nTl /2 . (44) 

The time evolution of the temperature is determined by Eq. (|38|), i.e. by the 
time dependence of /j 2 . 

The time-dependent coefficients /i p (t) may be expressed in terms of a 2 
owing to their definition Eq. ( |3~^ ) and the approximation / = 0(c) [1 + 
o 2 (t)5 2 (c 2 )]. We finally obtain: 



\ J dCl J dC2 j de ®(~ Cl2 ' e )l c i2 ' e|0(ci)<^(c 2 ) 

x {1 + a 2 [S 2 (c 2 ) + S2(c 2 2 )] + a\ S 2 (c 2 )S 2 (c 2 )} Z\(c? + c£) (45) 
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with the definition of A{c^ + g£) given above. Calculations performed up to 
the second order in terms of the dissipative parameter S yield 26 1: 

2 2 

t*2 = Y,Y, Akn6 ' ka % ( 4fi ) 



k=0 n=0 



where 



Aqo = 0; Aoi = 0; A02 = 

A w = w ; An = Au = 2500^° ' 

119 „ 4641 

Am = -wi : >42i = wi ; ^22 = wi 

400 640000 



with 



uj = 2\/2^2 1/10 r f^J Ci = 6.48562 .. . (48) 
wi = V2^2 1/5 r ^ C 2 = 9.28569 . . . (49) 

Similarly 

2 2 

M4 = ^]TW fe a 2 l (50) 



fc=0 n=0 



with 



Boo = 0; Bbi = 4^; £ 02 = 

56 „ 1806 567 

6x0 = ^0,0; B ^ = ^o; 612 = ( 51 ) 

77 149054 348424 

^ 20 - "lo"' 11 Bai - "T3750*' 15 622 _ "5500000 Wl 

Thus, Eqs. (||) and @, together with Eqs. @ and @ form a closed 
set to find the time evolution of the temperature and the coefficient a 2 . We 
want to stress an important difference for the time evolution of temperature 
for the case of the impact- velocity dependent restitution coefficient, as com- 
pared to that of a constant restitution coefficient. In the former case it is 
coupled to the time evolution of the coefficient 02, while in the latter case 
there is no such coupling since 02 = const. This coupling may lead in to a 
rather peculiar time-dependence of the temperature. 



Granular Gases with impact-velocity dependent restitution coefficient 



111 



Introducing the reduced temperature u(t) = T(t)/To we recast the set of 
equations (B8) and into the following form: 



5 2 7 2 

3 + 5° 2 + 500 &2 



-1 r 17/10 /' 5 119 1547 "I 



a 2 - r a u 1/2 n 2 (1 + a 2 ) + \r Q u 1/2 ii A = , (53) 

5 



where we introduce the characteristic time 



to 1 = jq S ■ T^O)- 1 (54) 



with 



q = 2 1/5 J T(21/10)Ci/8 = 5~ 2/5 Vtt^(3/5)/8 = 0.173318.. . (55) 

r ° s ik™" (56) 

gi = 2 1/10 {lo 1 /lu ) = 1.53445. . . (57) 

As shown below the characteristic time tq describes the time evolution of 
the temperature. To obtain these equations we use the expressions for fJ>2(t), 
B(t), and for the coefficients A n k- Note that the characteristic time To is 
S^ 1 ^> 1 times larger than the collision time ~ r c (0). 

We will find the solution to these equations as expansions in terms of the 
small dissipative parameter S (S'(t) = 6 ■ 2 1 / 10 u 1 / 10 (i)): 

u = uq + 8 ■ i*i + S 2 ■ u 2 + ■ ■ ■ (58) 
a 2 = a 2 o + 5 ■ a 2 i + 5 2 ■ a 22 + ■ ■ • (59) 

Substituting Eqs. @^H@ into Eqs. (H|||), one can solve these equa- 
tions perturbatively, for each order of S. The solution of the order of 0(1) 
reads for the coefficient a 2 (t) pq| : 

o ao (t) « a 20 (0) e - 4t /( 5rE (°» , (60) 

where te = §t c is the Enskog relaxation time, so that a 20 (t) vanishes for 
t ~ r . This refers to the relaxation of an initially non-Maxwellian veloc- 
ity distribution to the Maxwellian distribution. Note that the relaxation oc- 
curs within few collisions per particle, similarly to the relaxation of common 
molecular gases. 

We now assume that the initial distribution is Maxwellian, i.e., that 
O20 (0) = for t = 0. Then the deviation from the Maxwellian distribution 
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originates from the inelasticity of the interparticle collisions. For the case of 
fl 2o(0) = the solution of the order of 0(1) for the reduced temperature 
reads 

",•(" f I • -) , (61) 



To \ r 

which coincides with the time-dependence of the temperature obtained pre- 
viously using scaling arguments |23[ (up to a constant tq which may not be 
determined by scaling arguments). 

The solution for a%{t) in linear approximation with respect to S reads 

oa(i) = S ■ oat (t) = —wit)- 1 {Li [w(t)] - Li HO)]} (62) 

where 

w(t) = exp [(qoSy 1 (1 + t/r ) 1/6 ] . (63) 
and with the logarithmic Integral 

X 

u ^=Swf- (64) 



For t <C r the coefficient a^{t) ( p2[ ) reduces to 

a 2 (i)--^-/i(l-e- 4 */( 5TE (°») (65) 

where 

h = 2 1/10 (Bio - 5Ao) /16tt = (3/10)r(21/10)2 1/5 Ci = 0.415964 . (66) 



As it follows from Eq. (p5J) , after a transient time of the order of few collisions 
per particle, i.e. for te(0) < t -C To, Oa(t) saturates at the "steady-state" - 
value — /i<5 = —0.415964(5, i.e. it changes only slowly on the time-scale ~ 
t c (0). On the other hand, for t 3> To one obtains 

a 2 (t)~-6-h(t/ To )- 1/6 (67) 

so that 02(f) decays to zero on a time-scale ~ To, i.e. slowly in the collisional 
time-scale ~ t c (0) <C t . The velocity distribution thus tends asymptotically 
to the Maxwellian distribution. For this regime the first-order correction for 
the reduced temperature, u±(t), reads p(if : 

ux{t) = {^h + 2q 1 S \ (t/ro)- 11 / 6 = 3.26856 (t/r )- n / 6 , (68) 
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where we used the above results for the constants h and q\. From the last 
equation one can see how the coupling between the temperature evolution 
and the evolution of the velocity distribution influences the evolution of tem- 
perature. Indeed, if there were no such coupling, there would be no coupling 
term in Eq. (f>^), and thus no contribution from to the prefactor of Ui(t) 
in Eq. (|68|). This would noticeably change the time behavior of Ui(t). On the 
other hand, the leading term in the time dependence of temperature, Ua(t), 
is not affected by this kind of coupling. 

In Fig. H and Fig. |^ we show the time dependence of the coefficient d2(t) 
of the Sonine polynomial expansion and of the temperature of the granular 
gas. The analytical findings are compared with the numerical solution of 
the system ( |5^ , |53| ). As one can see from the figures the analytical theory 
reproduces fairly well the numerical solution for the case of small S. 




5 10 15 20 1 100 10000 

time time 




Time 



Fig. 2. Time dependence of the second coefficient of the Sonine polynomial ex- 
pansion et2(i). Time is given in units of the mean collisional time t c (0). (Left): 
a 2 x 1000 (solid lines) for S = 0.001,0.005,0.01,0.015 (top to bottom) together 
with the linear approximation (dashed lines); (Right): the same as (left) but for 
larger times; (Middle): -a 2 (t) over time (log-scale) for S = 0.03,0.01,0.003,0.001 
(top to bottom) together with the power-law asymptotics ~ t -1 ' 6 . 
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Fig. 3. Time-evolution of the reduced temperature, u(t) = T(t)/To. The time is 
given in units of mean collisional time r c (0). Solid line: numerical solution, short- 
dashed: Uo(t) = (1 + t/ro)~ B ^ 3 (zero-order theory), long-dashed: u(t) = uo(t) + 
8ui(t) (first-order theory). (Left): for S = 0.05,0.1 (top to bottom); (Right): 5 — 
0.15,0.25 (top to bottom); (Middle): the same as (Left) but log-scale and larger 
ranges. 



As it follows from Fig. || (where the time is given in collisional units), 
for small S the following scenario of evolution of the velocity distribution 
takes place for a force-free granular gas. The initial Maxwellian distribu- 
tion evolves to a non-Maxwellian distribution, with the discrepancy between 
these two characterized by the second coefficient of the Sonine polynomials 
expansion a-z. The deviation from the Maxwellian distribution (described by 
(12) quickly grows, until it saturates after a few collisions per particle at a 
"steady-state" value. At this instant the deviation from the Maxwellian dis- 
tribution is maximal, with the value 02 ~ —0.4(5 (Fig. 0a). This refers to 
the first "fast" stage of the evolution, which takes place on a mean-collision 
time-scale ~ r c (0). After this maximal deviation is reached, the second "slow" 
stage of the evolution starts. At this stage 0,2 decays to zero on the "slow" 
time scale r ~ 5 _1 t c (0) 3> tq(0), which corresponds to the time scale of 
the temperature evolution (Fig. 0b); the decay of the coefficient a% (t) in this 
regime occurs according to a power law ~ t^ 1 ^ 6 (Fig. 0c). Asymptotically 
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the Maxwellian distribution would be achieved, if the clustering process did 
not occur. 

Fig. H illustrates the significance of the first-order correction Ui(t) in the 
time-evolution of temperature. This becomes more important as the dissipa- 
tion parameter 6 grows (Figs. ||a,b). At large times the results of the first- 
order theory (with u\(t) included) practically coincide with the numerical 
results, while zero-order theory (without ui(t)) demonstrates noticeable de- 
viations (Fig. ||c). 

For larger values of <5 the linear theory breaks down. Unfortunately, the 
equations obtained for the second order approximation 0(8 2 ) are too com- 
plicated to be treated analytically. Hence, we studied them only numerically 
(see Fig. |^). As compared to the case of small 6, an additional intermediate 
regime in the time-evolution of the velocity distribution is observed. The first 
"fast" stage of evolution takes place, as before, on the time scale of few colli- 
sions per particle, where maximal deviation from the Maxwellian distribution 
is achieved (Fig. [|). For 8 > 0.15 these maximal values of 02 are positive. 
Then, on the second stage (intermediate regime), which continues 10 — 100 
collisions, a-i changes its sign and reaches a maximal negative deviation. Fi- 
nally, on the third, slow stage, 02 (t) relaxes to zero on the slow time-scale 
~ To, just as for small 6. In Fig. ^ we show the first stage of the time evo- 
lution of a,2(t) for systems with large 5. At a certain value of the dissipative 
parameter S the behavior changes qualitatively, i.e. the system then reveals 
another time scale as discussed above. 




5 10 15 20 

time 



Fig. 4. Time dependence of the second coefficient of the Sonine polynomial ex- 
pansion ct2(i) x 100. Time is given in units of mean collisional time t c (0). 8 = 
0.1, 0.11, 0.12, ... , 0.20 (bottom to top). 

Figure |H shows the numerical solution of Eqs. ( p2| ) and ( p3| ) for the second 
Sonine coefficient 02 (t) as a function of time. One can clearly distinguish the 
different stages of evolution of the velocity distribution function. A more 
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detailed investigation of the evolution of the distribution function for larger 
dissipation is subject of present research [p6[. 




The interesting property of the granular gases in the regime of homoge- 
neous cooling is the overpopulation of the high-velocity tails in the velocity 
distribution which has been shown for granular gases consisting of par- 
ticles which interact via a constant restitution coefficient, e = const. How 
does the velocity dependence of the restitution coefficient as it appears for 
viscoelastic spheres influence this effect? We observe, that for the case of 
e = e(wimp) the functional form (i.e. the exponential overpopulation ||) per- 
sists, but it decreases with time on the "slow" time-scale ~ to. Namely we 
obtain for the velocity distribution for c » 1 Q : 



/(c, t) ~ exp 



b ( t' 



1/61 



(69) 



where b = (16q / 5) 1 = 2.25978 . . . , which holds for t > r c (0). Again 

we see that the distribution tends asymptotically to the Maxwellian distri- 
bution, since the overpopulation vanishes as t — > oo. 



Granular Gases with impact-velocity dependent restitution coefficient 



117 



Using the temperature and the velocity distribution of a granular gas as 
were derived in this section, one can calculate the kinetic coefficients. In the 
next section we consider the simplest one - the self-diffusion coefficient. 

4 Self-diffusion in granular gases of viscoelastic 
particles 

In the simplest case diffusion of particles occurs when there are density gra- 
dients in the system. The diffusion coefficient D relates the flux of particles 
J to the density gradient Vrt according to a linear relation, provided the 
gradients are not too large: 

J = -D\7n. (70) 

The coefficient D also describes the statistical average of the migration of a 
single particle. For equilibrium 3D-systems the mean-square displacement of 
a particle reads 

((Ar(t)f) =6Dt, (71) 

\ / eq 

where (• • • ) denotes the equilibrium ensemble averaging. For nonequilibrium 
systems, such as granular gases, one should consider the time-dependent dif- 
fusion coefficient D(t) and the corresponding generalization of Eq. (|7l|): 

((Ar(t)) 2 ) = 6 ^D^dt', (72) 

where (• • • ) denotes averaging over the nonequilibrium ensemble. If the mi- 
gration of a particle occurs in a uniform system composed of particles of the 
same kind, this process is called "self-diffusion" . Correspondingly, the kinetic 
coefficient D is called self-diffusion coefficient. 

To find the mean-square displacement, one writes 

((Ar(t)f) = (J* v{t')dt> £ «(«")&") (73) 
and encounters then with the velocity autocorrelation function 

K v (t',t) = {v(t')v(t")) 

which should be evaluated in order to obtain the mean-square displacement 
and the self-diffusion coefficient. 

To calculate K v (t', t) we use the approximation of uncorrelated successive 
binary collision, which is valid for moderately dense systems, and an approach 
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based on the formalism of the pseudo-Liouville operator C |27j] . The pseudo- 
Liouville operator is defined as 

^ = £*i-^T+E% ( 74 ) 



7 

1<J 



The first sum in ( [74] ) refers to the free streaming of the particles (the ideal 
part) while the second sum refers to the particle interactions which are de- 
scribed by the binary collision operators p3] 

f ij =a 2 jdee( c v ir e)\v ij -e\8{r ij -ae){b^-l) , (75) 
where 0(x) is the Heaviside function. The operator 6|- is defined as 

%f ( r i> T h V h V 3---) = f ( r n ' ■ ') > ( 76 ) 

where / is some function of the dynamical variables and v* and v* are the 
postcollisional velocities from Eq. The pseudo-Liouville operator gives 

the time derivative of any dynamical variable B (e.g. p4|| ): 

j t B ({n, Vl } ,t) = iCB ({n, v { } , t) . (77) 
Therefore, the time evolution of B reads (t > t') 

B ({ri, Vi },t) = e iC ^B ({r,, «,},*') . (78) 
With Eq. (fsl) the time-correlation function reads 

(B(t')B(t)) = f dr P (t')B(t')e lC ( t - t ">B(t') , (79) 



where J dr denotes integration over all degrees of freedom and p(t') depends 
on temperature T, density n, etc., which change on a time-scale t > r c . 
Now we assume that 

(i) the coordinate part and the velocity part of the distribution function p(t) 
factorize, and 

(ii) the molecular chaos hypothesis is valid. 

This suggests the following form of the distribution function: 

p(t) = p(n, ... ,r N )- f(v u t) . . . f(v N , t) . (80) 

In accordance with the molecular chaos assumption the sequence of the suc- 
cessive collisions occurs without correlations. If the variable B does not de- 
pend on the positions of the particles, its time-correlation function decays 
exponentially p9|: 
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(B(t')B(t)) = (B 2 ) t , e~ l*-*'IA B (*') (t > t') . (81) 

where (•••)*' denotes the averaging with the distribution function taken at 
time t' . The relaxation time tb is inverse to the initial slope of the autocorre- 
lation function , as it may be found from the time derivative of (B(t')B(t)) 
taken at t = t'. Equations (|7^) and (^lj) then yield 

-T B \t>) = Jdr P (t')BiCB/(B\ = (82) 

The relaxation time T^it') depends on time via the distribution function 
p(t') and varies on the time-scale t S> t c . 

Let B(t) be the velocity of some particle, say Vi(t). Then with 3T(t) = 
(v 2 ) v Eqs. © and © (with Eqs. @ and ©) read JIo| 

(171(0 ' »i(*)> = ST^Oe-l*-*'!/^*') (83) 



«i • T 12 v x 

r-t(t') = (N-l)±- ^1. (84) 

{vi 

To obtain Eq. j8^ ) we take into account that CqV\ = 0, TyVi — (for i ^ 1) 
and the identity of the particles. 

Straightforward calculation yields for the case of a constant restitution 
coefficient: 

r* l (t) = ^ln<? 2 g2(<?)V^nt)= e -^T E \t), (85) 

where teU) — § T c (t) is the Enskog relaxation time [f24[ . Note that according 
to Eq. (pq), t„ = j^te > tb, i.e., the velocity correlation time for inelastic 
collisions exceeds that of elastic collisions. This follows from partial suppres- 
sion of the backscattering of particles due to inelastic losses in their normal 
relative motion, which, thus, leads to more stretched particle trajectories, as 
compared to the elastic case. 

Similar (although somewhat more complicated) computations may be per- 
formed for the system of viscoelastic particles yielding 



T« 1 (t)=T E -\t) 



l + ^a 2 (t)-iq 5u^ w (t) 
16 



(86) 



where qo — 0.173318 has been already introduced and a%{t), u(t) are the 
same as defined above. To obtain Eq. (^6|) we neglect terms of the order of 
0(4), 0{8 2 ) and 0(a 2 S). 



120 N. V. Brilliantov, T. Poschel 



Using the velocity correlation function one writes 



((Ar(t)) 2 ^ = 2 ( dt'3T(t') I dt" e - | *"-*' l/T « (t ' ) . (87) 



On the short-time scale t ~ r c , T(t') and r v (t') may be considered as con- 
stants. Integrating in (j87j) over t" and equating with ([72|) yields for £ ^> r c ~ 
t„ the time-dependent self-diffusion coefficient 

D(t)=T{t)T v (t). (88) 

Substituting the dependencies for u(t) = T(t)/To and a%{t) as functions 
of time, which has been derived in the previous section, we obtain the time 
dependence of the coefficient of self-diffusion D(t). For t 3> tq this may be 
given in an explicit form: 



D \toJ V 400 J \t 

where the constants qo, q\ and h are given above. Hence, the prefactor in 
the term proportional to S reads (4qo + qi + = 2.24956, and D is the 

initial Enskog value of the self-diffusion coefficient 

= \^ng 2 {a)aV 2 . (90) 

Correspondingly, the mean-square displacement reads asymptotically for t ^> 
To: 

({Ar{t)f\ ~t x ' 6 + b8\ogt + ... , (91) 

where 6 is some constant. This dependence holds true for times 

t^S- 1 <i< r c (0) <T n / 5 . (92) 



The first inequality in Eq. (|92|) follows from the condition tq <C t, while 
the second one follows from the condition T c (t) <C To, which means that 
temperature changes are slow on the collisional time-scale. For the constant 
restitution coefficient one obtains 

T(t)/T = [l + 7 ot/r c (0)r 2 , (93) 

where 70 = (l - e 2 ) /6 [|§. Thus, using Eqs. ( |85| ) and ( |S8| ) one obtains for 
the mean-square displacement in this case 

((Ar(t)) 2 )~logt. (94) 
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As it follows from Eqs. (^Tj) and ( f)4| ) the impact- velocity dependent restitu- 
tion coefficient, Eq. (10), leads to a significant change of the long-time be- 
havior of the mean-square displacement of particles in the laboratory-time. 
Compared to its logarithmically weak dependence for the constant restitution 
coefficient, the impact-velocity dependence of the restitution coefficient gives 
rise to a considerably faster increase of this quantity with time, according to 
a power law. 

One can also compare the dynamics of the system in its inherent-time 
scale. First we consider the average cumulative number of collisions per par- 
ticle Af(t) as an inherent measure for time (e.g. 0,f§). It may be found by 
integrating dftf — T c (t) _1 <ii ||. For a constant restitution coefficient e one 
obtains Af(t) ~ logi, while for the impact- velocity dependent e (t>i mp ) one has 
Af(t) ~ t 1 / 6 . Therefore, the temperature and the mean-square displacement 
behave in these cases as 



e = const 


e = e (w imp ) 


T(A0 ~ e- 2 ^ 2 )* 


T(Af ) ~ N- 10 







If the number of collisions per particle M{t) would be the relevant quantity 
specifying the stage of the granular gas evolution, one would conjecture that 
the dynamical behavior of a granular gas with a constant e and velocity- 
dependent e are identical, provided an AT-based time-scale is used. Whereas 
in equilibrium systems the number of collisions is certainly an appropriate 
measure of time, in nonequilibrium systems this value has to be treated with 
more care. As a trivial example may serve a particle bouncing back and 
forth between two walls, each time it hits a wall it loses part of its energy: 
If one describes this system using a A/"-based time, one would come to the 
conclusion that the system conserves its energy, which is certainly not the 
proper description of the system. According to our understanding, therefore, 
the number of collision is not an appropriate time scale to describe physical 
reality. Sometimes, it may be even misleading. 

Indeed, as it was shown in Ref. ||, the value of J\f c , corresponding to a 
crossover from the linear regime of evolution (which refers to the homoge- 
neous cooling state) to the nonlinear regime (when clustering starts) may 
differ by orders of magnitude, depending on the restitution coefficient and on 
the density of the granular gas. Therefore, to analyze the behavior of a granu- 
lar gas, one can try an alternative inherent time-scale, T _1 = T(t)/To which 
is based on the gas temperature. Given two systems of granular particles 
at the same density and the same initial temperature To, consisting of par- 
ticles colliding with constant and velocity-dependent restitution coefficient, 
respectively, the time T allows to compare directly their evolution. A strong 
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argument to use a temperature-based time has been given by Goldhirsch and 
Zanetti || who have shown that as the temperature decays, the evolution of 
the system changes from a linear regime to a nonlinear one. Recent numerical 
results of Ref. ^ also support our assumption: It was shown that while Af c 
differs by more than a factor of three for two different systems, the values 
of T c , (defined, as T c = T(J\f c )/T ) are very close [pj. These arguments show 
that one could consider T as a relevant time-scale to analyze the granular 
gas evolution. 

With the temperature decay T(Af)/To ~ e" 270 ^ for a constant restitu- 
tion coefficient and T(Af)/To ~ JV~ 10 for the impact-velocity dependent one, 
we obtain the following dependencies: 



e = const 


e = e (fimp) 


T ~ ± 

1 T 




T ~ -1 


((Ar(T)f) 


~ logT 


UAr(T)fS ~ T l l w 



This shows that in the temperature-based time-scale, in which the cool- 
ing of both systems is synchronized, the mean-square displacement grows 
logarithmically slow for the case of constant restitution coefficient and much 
faster, as a power law, for the system of viscoelastic particles with e = e{n lmp ). 
Thus, we conclude that clustering may be retarded for the latter system. 



5 Conclusion 



In conclusion, we considered kinetic properties of granular gases composed 
of viscoelastic particles, which implies the impact-velocity dependence of the 
restitution coefficient. We found that such dependence gives rise to some 
new effects in granular gas dynamics: (i) complicated, non-monotonous time- 
dependence of the coefficient ai of the Sonine polynomial expansion, which 
describes the deviation of the velocity distribution from the Maxwellian and 
(ii) enhanced spreading of particles, which depends on time as a power law, 
compared to a logarithmically weak dependence for the systems with a con- 
stant e. 

The Table below compares the properties of granular gases consisting of 
particles interacting via a constant coefficient of restitution e = const and 
consisting of viscoelastic particles where the collisions are described using an 
impact velocity dependent restitution coefficient e = e(i>imp): 
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e = const 


e = e («imp) 


e is a model parameter 


e = i-c 1 ^ 2 /^ 5 p + ... 

Ci = 1.15396, C 2 = fCV 2 , . . . 
K = k(Y, v, m, R) 
A = A( m , m ,Y,y 

all quantities are defined via parameters 
of the particle material Y, v, i] 1 / 2 
and their mass and radius. 


Small 

1 — e 2 - does not depend 
on the state of the system 


parameter 

5 = An 2 ' b Tl /W - depends 
on the initial temperature Tq. 


Temperature 

t = t (i + t/r^y 2 \t = T (i + t/T y 5/3 


Velocity 

/(«.*) = ^ij/» 

f(c,t) = ^c){l + ^ =1 a p S p (c 2 )} 

ci2 = const. 


distribution 

/»=0(c){i+Er=i«pW^( c2 )} 

a2 = ci2(t) - is a (complicated) 
function of time. 


Self- 

((Zk(i)) 2 )~logt 


diffusion 

((Ar(t)) 2 ) ~ tVe 
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